Matlab Codes

Chapter 2

Matlab code 2.6: Matlab file “Figure 2-11.m”

%--------------------------------------------------------------------

% This code can be used to generate Figure 2.11

%--------------------------------------------------------------------

 

clear all;

close all;

%% the theoretical curve of precession target micro-Doppler characteristic when radar platform is vibrating

c = 3e8;

j = sqrt(-1);

fc = 10e9; % carrier frequency of transmitted signal

v = 0; % translational velocity of target

cord = 1000*[300 100 500]; % coordinates of local coordinate system's origin in the radar coordinate system

colo = [0 0 0]; % coordinates of radar in the radar coordinate system

R0 = cord-colo;

n = R0/sqrt(sum(R0.^2)); % unit vector of Line-of-Sight(LOS)

p1 = [0 0 0.5]; % scatterer on the apex of cone

p2 = [0.5 0 -0.5]; % scatterer on the bottom of cone

tgt = [p1;p2];

alpha = atan(cord(2)/cord(1)); % azimuth angle

beita = atan(cord(3)/(sqrt(cord(1)^2+cord(2)^2))); % elevation angle

alpha0 = pi/10; % azimuth angle of radar platform vibrating axis

beita0 = pi/3; % elevation angle of radar platform vibrating axis

dr = 0.01; % vibration amplitude of radar platform

fr = 8; % vibration frequency of radar platform

ae = pi*[1/3 1/4 1/5]; % Euler angle

ws = [19.2382 -11.1072 22.2144]; % angular velocity of spinning in the local coordinate system

omegas = sqrt(sum(ws.^2));

% omegas = 10*pi;

Ts = 0.2; % spinning period

wc = [7.6756 -2.3930 9.6577]; % angular velocity of coning in the local coordinate system

omegac = sqrt(sum(wc.^2));

% omegac = 4*pi;

Tc = 0.5; % coning period

t = 2;% radar illumimated time

ri1 = [cos(ae(3)) sin(ae(3)) 0;-sin(ae(3)) cos(ae(3)) 0;0 0 1];

ri2 = [1 0 0;0 cos(ae(2)) sin(ae(2));0 -sin(ae(2)) cos(ae(2))];

ri3 = [cos(ae(1)) sin(ae(1)) 0;-sin(ae(1)) cos(ae(1)) 0;0 0 1];

ri = ri1*ri2*ri3; % initial rotation matrix

ws = ri*ws'; % angular velocity of spinning in the reference coordinate system

% omega = sqrt(sum(w.^2));

wse = ws/omegas; % unit vector of spinning

wsr = [0 -wse(3) wse(2);wse(3) 0 -wse(1);-wse(2) wse(1) 0]; % skew symmetric matrix of spinning

wc = ri*wc'; % angular velocity of coning in the reference coordinate system

wce = wc/omegac; % unit vector of coning

wcr = [0 -wce(3) wce(2);wce(3) 0 -wce(1);-wce(2) wce(1) 0]; % skew symmetric matrix of coning

r0a = tgt(1,:)-colo;

r0b = tgt(2,:)-colo;

r0ar = ri*r0a'; % spinning scatterer in the reference coordinate system

r0br = ri*r0b'; % coning scatterer in the reference coordinate system

a1 = (2*fc/c)*r0br'*(omegac*wcr^2+omegac*wcr^2*wsr^2)'*n';

a2 = (2*fc/c)*r0br'*(omegac*wcr+omegac*wcr*wsr^2)'*n';

a3 = (2*fc/c)*r0br'*(omegac*wcr^2*wsr+omegas*wcr*wsr^2)'*n';

a4 = (2*fc/c)*r0br'*(omegas*wcr*wsr-omegac*wcr^2*wsr^2)'*n';

a5 = (2*fc/c)*r0br'*(omegac*wcr*wsr-omegas*wcr^2*wsr^2)'*n';

a6 = (2*fc/c)*r0br'*(-omegac*wcr*wsr^2-omegas*wcr^2*wsr)'*n';

a7 = (2*fc/c)*r0br'*(omegas*wsr^2+omegas*wcr^2*wsr^2)'*n';

a8 = (2*fc/c)*r0br'*(omegas*wsr+omegas*wcr^2*wsr)'*n';

prf = 4000; % pulse repetition frequency

pri = 1/prf; % pulse repetition interval

dt = 0:pri:t-pri; % time sampling interval

m = length(dt);

fmd1 = zeros(2,m); % theoretical micro-Doppler frequency

fmd2 = zeros(2,m);

for i = 1:m

    fmd1(1,i) = a1*sin(omegac*dt(i))+a2*cos(omegac*dt(i))+a3*sin(omegas*dt(i))*sin(omegac*dt(i))...

               +a4*sin(omegac*dt(i))*cos(omegas*dt(i))+a5*cos(omegac*dt(i))*sin(omegas*dt(i))...

               +a6*cos(omegac*dt(i))*cos(omegas*dt(i))+a7*sin(omegas*dt(i))+a8*cos(omegas*dt(i))...

               -(4*pi*fc*fr*dr/c)*(cos(alpha-alpha0)*cos(beita)*cos(beita0)+sin(beita)*sin(beita0))*cos(2*pi*fr*dt(i));

    fmd1(2,i) = (2*fc*omegac/c).*((wcr^2.*sin(omegac*dt(i))+wcr.*cos(omegac*dt(i)))*ri*r0ar)'*n'...

               -(4*pi*fc*fr*dr/c)*(cos(alpha-alpha0)*cos(beita)*cos(beita0)+sin(beita)*sin(beita0))*cos(2*pi*fr*dt(i));

    fmd2(1,i) = a1*sin(omegac*dt(i))+a2*cos(omegac*dt(i))+a3*sin(omegas*dt(i))*sin(omegac*dt(i))...

               +a4*sin(omegac*dt(i))*cos(omegas*dt(i))+a5*cos(omegac*dt(i))*sin(omegas*dt(i))...

               +a6*cos(omegac*dt(i))*cos(omegas*dt(i))+a7*sin(omegas*dt(i))+a8*cos(omegas*dt(i));

    fmd2(2,i) = (2*fc*omegac/c).*((wcr^2.*sin(omegac*dt(i))+wcr.*cos(omegac*dt(i)))*ri*r0ar)'*n';

end

figure(1)

plot(dt,fmd1,'r') % when radar platform is vibrating

hold on

plot(dt,fmd2,'b') % when radar platform is stationary

xlabel('Time (s)')

ylabel('Frequency (Hz)')

axis([0,2,-1200,1200])

 

%% the time-frequency analysis result of precession target micro-Doppler characteristic when radar platform is vibrating

r = zeros(length(tgt(:,1)),length(dt)); % distance between the scatterers and radar

for i = 1:m

    for n = 1:length(tgt(:,1))

        if n == 1

           r(n,i) = sqrt((R0'+v*t+((eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))*r0ar)...

                        -(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)]))'...

                        *(R0'+v*t+((eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))*r0ar)...

                        -(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)])));

        else

           r(n,i) = sqrt((R0'+v*t+(eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))...

                        *(eye(3)+wsr*sin(omegas*dt(i))+wsr^2*(1-cos(omegas*dt(i))))*r0br...

                        -(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)]))'...

                        *(R0'+v*t+(eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))...

                        *(eye(3)+wsr*sin(omegas*dt(i))+wsr^2*(1-cos(omegas*dt(i))))*r0br...

                        -(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)])));

        end

    end

end

s = zeros(length(dt),length(tgt(:,1))); % echo signal matrix

for i = 1:length(tgt(:,1))

    s(:,i) = exp(j*2*pi*fc*2*r(i,:)'/c);

end

st = sum(s,2);

N = 200; % number of Gabor coefficients in time

Q = 50; % degree of oversampling

tfr = tfrgabor(st.',N,Q); % Gabor transform

tfr_r = fftshift(tfr,1);

figure(2)

contour(linspace(min(dt),max(dt),length(tfr_r(1,:))),linspace(-1000,1000,length(tfr_r(:,1))),tfr_r,70)

xlabel('Time (s)')

ylabel('Frequency (Hz)')

axis([0,2,-800,800])